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The characteristic approach to numerical relativity is a useful tool in evolving gravitational sys- 
tems. In the past this has been implemented using two patches of stereographic angular coordinates. 
In other applications, a six-patch angular coordinate system has proved effective. Here we investi- 
gate the use of a six-patch system in characteristic numerical relativity, by comparing an existing 
two-patch implementation (using second-order finite differencing throughout) with a new six-patch 
implementation (using either second- or fourth-order finite differencing for the angular derivatives). 
We compare these different codes by monitoring the Einstein constraint equations, numerically eval- 
uated independently from the evolution. We find that, compared to the (second-order) two-patch 
code at equivalent resolutions, the errors of the second-order six-patch code are smaller by a factor 
of about 2, and the errors of the fourth-order six-patch code are smaller by a factor of nearly 50. 

PACS numbers: 04.25.Dm 



I. INTRODUCTION 

The characteristic, or null-cone, approach to numerical relativity is based on the Bondi-Sachs metric and has 

been successfully implemented in the PITT code H S 0, H 0, 1 (and also in other codes H IH IH llllll El 113). 
The PITT code has the important property that all tests have shown it to be long-term stable, for example in evolving 
a single black hole spacetime |l6j . However, there are problems with the computation of gravitational radiation: it 
has been successfully computed in test cases 0, Hi and also in scattering problems ^3 , but it has not been possible to 

\ compute gravitational radiation emitted in astrophysically interesting scenarios, such as a star in close orbit around 

. a black hole 0,111. 

\ This paper is a contribution towards the longer term goal of producing a characteristic code that can reliably 
compute gravitational radiation in situations of astrophysical interest, either within a stand-alone code 0,^3, or 
within the context of Cauchy-characteristic extraction or Cauchy-characteristic matching [l9l | . The first step in this 
process was to find a class of exact solutions to the linearized Einstein equations |2Cll |. These solutions are written 
in terms of the Bondi-Sachs metric and continuously emit gravitational radiation. Thus they provide an appropriate 

[ testbed for validating further developments of the characteristic code. 

In this paper, we investigate a strategy for improving the code. In another context in numerical relativity |2ll |. 

\ good results have been obtained by coordinatizing the sphere by means of six angular patches, rather than the two 
patch stereographic coordinates used in the PITT code. One advantage of the six patch system is that inter-patch 
interpolation is particularly simple because it is one-dimensional (since ghost zone points lie on grid lines), leading 
to an expectation of reduced noise at the patch interfaces. Further, this means that it is easy to go from second- to 
fourth-order accurate finite differencing of angular derivatives, since the interpolation order must also be increased 
and this is straightforward in the one-dimensional case. Using the exact solutions of ref. [23| as a testbed, we have 
computed and compared the errors obtained in (a) the stereographic code, (b) the second-order six patch code, and (c) 
the fourth-order six patch code. Note that, in the fourth-order case, the radial and time derivatives are second-order 
accurate, so that overall the code is second-order - unfortunately, changing these derivatives to fourth-order would be 
a major undertaking, because they arc not confined to a single subroutine. We found that, for equivalent resolutions, 
the errors of (b) are a little smaller than those of (a) . and that the errors of (c) are much smaller than those of (a) by 
a factor of nearly 50. 

As in the ADM formalism, four of the ten characteristic Einstein equations are not used in the evolution but 
constitute constraints. We have constructed and validated code that evaluates the constraints, and this can be used 
as a tool to monitor the reliability of a computational evolution. 

The numerical computations presented in this paper were all performed within the Cactus computational 
toolkit |22l | (http://www.cactuscode.org), using the Carpet driver (http:V7wwwTcarpetcode.org) to support 
the multiple-patch computations. The computer algebra results were obtained using Maple. 

The plan of the paper is as follows. Sec. ^ summarizes background material that will be used later. Sec. Illll 
describes our implementation of the six patch angular coordinate system. Sec. lIVI describes the constraint evaluation. 
Computational results are presented in Sec. |3 and are then discussed in the Conclusion, Sec. IVII 
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II. BACKGROUND MATERIAL 
A. The Bondi-Sachs metric 

The formalism for the numerical evolution of Einstein's equations, in null cone coordinates, is well known P, 0, S 
IMIEI^' For the sake of completeness, we give a summary of those aspects of the formalism that will be used here. 
We start with coordinates based upon a family of outgoing null hypersurfaces. We let u label these hypersurfaces, 
{A = 2,3), label the null rays and r be a surface area coordinate. In the resulting = (u, r, x^) coordinates, the 
metric takes the Bondi-Sachs form 0, Q 

ds^ = - (e2^(l + W,r) ~ r^hABU'^U") du^ 
-2e^'^dudr ~ 2r^hABU^dudx^ + r^hAsdx^dx^, (1) 

where h^^hBc = ^'^'^ det{hAB) = det{qAB), with qab a metric representing a unit 2-sphere embedded in flat 
Euclidean 3-space; Wc is a normalised variable used in the code, related to the usual Bondi-Sachs variable V by 
V = r + Wct"^ . As discussed in more detail below, we represent qab by means of a complex dyad qA- Then, for an 
arbitrary Bondi-Sachs metric, Hab can then be represented by its dyad component 

J = hABq^q^'l^, (2) 

with the spherically symmetric case characterized by J = 0. We also introduce the spin-weighted field 

U = f/^QA, (3) 

as well as the (complex differential) cth operators 9 and 3 j25l |. 

Einstein's equations Rap = 87r(TQ^ — \gapT) are classified as: hypersurface equations - Rn^q^RiAjh"^^ Rab - 
forming a hierarchical set for P,U and Wc] evolution equation q^q^RAB for J; and constraints Roa- An evolution 
problem is normally formulated in the region of spacetime between a timelike or null worldtube and future null infinity, 
with (free) initial data J given on u = 0, and with boundary data for (3, U, Wc, J satisfying the constraints given on 
the inner worldtube. 

B. The spin- weighted formalism and the 9 operator 

A complex dyad is written 

qA = {r2e''^',r3^^') (4) 

where rA,4'A are real quantities (but in general they are not vectors). The real and imaginary parts of qA are unit 
vectors that are orthogonal to each other, and qA represents the metric. Thus 

q^qA = 0, q'^qA = 2, qAB = ^{qAQB + qAqB)- (5) 

It is straightforward to substitute the a 2-metric into Eq. ((Sj to find ta and {4>^ — 4'2)- Thus qA is not unique, up to 
a unitary factor: if qA represents a given 2-metric, then so does = e*"gA. Thus, considerations of simplicity are 
used in deciding the precise form of dyad to represent a particular 2-metric. For example, the dyads commonly used 
to represent some unit sphere metrics, namely spherical polars and stercographic, are 

ds' = de'+sinH'd<t>' : qA = (l,zsin0); ds^ = ^C^^g^ +'^P^) . ^ 2 

Having defined a dyad, we may construct complex quantities representing all manner of tensorial objects, for 
example Xi = TAq^, X2 ~ T^^qAqB-, X3 = T^^ qAqB<f^ ■ Each object has no free indices, and has associated with 
it a spin-weight s defined as the number of q factors less the number of q factors in its definition. For example, 
s(Xi) = 1, 5(^2) = 0, 5(^3) = —3, and, in general, s{X) = —s{X). We define derivative operators S and 8 acting on 
a quantity V with spin-weight s 

W = q^dAV + sTV, W = q'^dAV - sfV (7) 
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where the spin- weights of dV and dV arc s + 1 and s — 1 , respectively, and where 

r = -ig-^q^V^te. (8) 

In the case of spherical polars, F = — cot 6, and for stcrcographic coordinates T = q + ip. 
The spin-weights of the quantities used in the Bondi-Sachs metric are 

s{Wc) = s{P) = 0, s(J) = 2, s{J) = -2, .s{U) = 1, s{U) = -1. (9) 

We will be using spin- weighted spherical harmonics |2^ |23| using the formalism described in ■ K will prove 
convenient to use gZem rather than the usual sYem (the suffix s denotes the spin-weight) as basis functions, where 

sZe,n = ^ [sYtm + i-iy"sYi-„,] for m > 

sZirn = --J= " sYl-m\ for m < 

s^^O = 5X^0, (10) 

and note that [T^ uses the notation sRim rather than the sZim used here; we use a different notation to avoid any 
confusion with the Ricci tensor. In the case s = 0, the s will be omitted, i.e. Zim = oZtm- Note that the effect of the 
3 operator acting on Zem is 



dZe^ = ^£{e + l)iZ,,n, = - + + 2) 2^^™- (11) 



C. Linearized solutions 

A class of solutions, in Bondi-Sachs form, to the linearized Einstein equations in vacuum was presented in rcf. |2fll |. 
and we use these solutions to test the accuracy of the numerical evolutions described later. More specifically, the 
solutions to be used are those given in Sec. 4.3 of ref. [23| for the case of a dynamic spacetime on a Minkowski 
background. We write 



J = ^{e - !)£{£ + l)i£ + 2) 2Zirn^{Ji{r)e'''"), [/ = ^£(7+1) iZ,™3?([/,(r)e*''"), 

p = Ze^MPee'"!, = Ze^,^{Wci{r)e'n, (12) 

where (r), Ui{r), (3i, Wd^r) are in general complex, and taking the real part leads to cos(z^it) and sin(i^M) terms. The 
quantities f3 and Wc are real; while J and U are complex due to the terms d^Zim and dZim, representing different 
terms in the angular part of the metric. We require a solution that is well-behaved at future null infinity, and is 
well-defined for r > 2, at which surface we set the inner boundary. Wc find in the case £ = 2 



Mr) 
U2{r) 



24/3o + 3wCi - iiy^C2 Ci C2 



36 4r 12r3 
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36 r 2r2 4r4 

24w/3o - 3i/2Ci + i^-'Ca 3iiyCi - 6/3o - iiy^C2 ^"^€2 ivC2 Q 



WMr) = ^ + + ^ + (13) 

with the (complex) constants (3o, Ci and C2 freely specifiable. 
We find in the case £ = 3 

[32 = l3o 

GOPo + iii^Ci + i'^C2 Ci WC2 C2 



Mr) = 
Mr) = 



180 lOr 6r3 4r4 

-GOiiyPo + Siy^Ci - w^C2 2/3o Ci 2v^C2 ^wC2 C2 
180 '^~'^2^ 3^ ^ 4r4 ^ ~^ 



, , miv[3Q~iv^Ci + w'^C2 iuC'i - 2[3o + iy*C2 i2u^C2 ^iy^C2 5iyC2 3C2 
W.,ir) = + + __ + _. (14) 
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The emitted gravitational radiation, that is the news N, takes a simple form in the linearized limit when the metric 
satisfies Eq. (O 

N^n (^e'"" ^hm (^^^^±^J, - Y^^JLr^ + e'^^'Pe^ ^/{e ^ !)£{£+ 1)(£ + 2) ^Zi^. (15) 



For the cases I =2 and 3, 



£ = 2 : N = ^\ '-^e^'A 2^2™; ^ = 3 : TV = 3? f-^^e'"^' ] 2^3™- (16) 



D. Six- patch Angular Coordinates 

Given Cartesian coordinates {x, y, z), we define "inflated-cube" angular coordinates on each 2-sphere of constant u 
and r, 

/.( = rotation angle about the x axis = arctan(?//z) 

V = rotation angle about the y axis = arctan(x/z) (17) 
(j) = rotation angle about the z axis = arctan(?//.T) 

where all the arctangents are 4-quadrant based on the signs of x, y, and z. We then introduce 6 coordinate patches 
covering neighborhoods of the ±z, ±x, and ±y axes, with the angular coordinates = (p, a) in each patch defined 
as follows (see note p9|'): 

{{v, /i) in ±z patches 
{v, (j)) in ±x patches (18) 
(/i, 4>) in zty patches 

Notice that each patch's x^ coordinates are nonsingular throughout a neighborhood of the patch, and that the 
union of all the patches covers without coordinate singularities. The name "inflated-cube" comes from another 
way to visualize these patches and coordinates: Imagine an xyz cube with xyz grid lines painted on its face. Now 
imagine the cube to be flexible, and inflate it like a balloon, so it becomes spherical in shape. The resulting coordinate 
lines will closely resemble those for (/i, i^, 0) coordinates. 

We introduce ghost zones in the usual manner along the angular boundaries of each patch, and we refer to the non- 
ghost-zone part of each patch's numerical grid as the "nominal" grid. We size the patches so they overlap slightly, 
with each ghost-zone grid point lying in the nominal grid of some other patch. Figure ^ shows an example of a 
six-patch system of this type. 

We couple the patches together by interpolating the field variables from neighboring patches to each ghost-zone 
grid point. Notice that with the definition H18|l . the angular coordinate x^ perpendicular to an interpatch boundary 
is always common to both the adjacent patches. This means that the interpatch interpolation need only be done in 
1 dimension, parallel to the interpatch boundary. 

Each patch uses its own x^ = (p, a) dyad for the spin- weighted quantities; when interpolating data from one patch 
to another we transform the data (after interpolation) as described in section UlI Al 

III. IMPLEMENTATION OF SIX PATCH ANGULAR COORDINATES 
A. Spin- weighted formalism 

The unit sphere metric qab in each patch is written as 

ds^ = (1 - sin^ p sin^ a)^^ ^ cos^ cr dp'^ + cos^ p da"^ - ^ sin(2p) sin(2cr) dpda^ , (19) 

with respect to coordinates (p, cr) with range (— 7r/4, -l-7r/4) x (— 7r/4, -t-7r/4). A (simple) dyad representing Eq. l(T^ 
is 

/ (6'c + »6>,) cos g {9c -i9s) cos p \ ^ ( 9s + iOc „ 9^ - i9c \ , 
i^pi ' 1^0201 j' 1 =(,20=^«^^'2^^^^^^j' (20) 



FIG. 1: This figure shows a six-patch grid covering 5^ at an angular resolution of Ax'* — 5°. Each patch's nominal grid is 
shown with solid lines; the central patch's ghost-zone points are shown with solid dots. The ghost zone has a width of 2 grid 
points (suitable for 5-point angular finite difference molecules). 



where 

' 1 — sin p sin a / 1 + ^^^^ P 



2 ' V 2 

The factor F, defined in Eq. ||SJ) and needed for the evaluation of 9, works out to be 

cos^ p cos^ cr(sin p + sin cr) + (cos^ p — cos^ a) (sin a — sin p) 



(21) 



r = 

Wc cos a cos p 

pnQ^ n rriQ^ rr{ Qin n — Qin rr\ 4- ( pnQ'^ rr — PPiS"^ p) (sin p -f- siu (j) ^ ^ 



46's COS (T COS p 

In addition, we need to specify how spin-weighted quantities transform at inter-patch boundaries. Suppose that we 
have two patches "Old" and "New", with quantities in the two patches being denoted by means of suffices (o) and 
(AT), respectively. Define the Jacobian Jg from the old patch to the new patch 

Then the dyad q^-j has components in the new coordinates 

^to)[N] ^ 1(0)^3 (24) 
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where the notation is that the (o) indicates that we are referring to the dyad that generates the metric in the old 
patch, and the [^v] indicates that the components are given with respect to the new coordinates. Any two complex 
dyads are related by means of a rotation exp(i7). Writing 

= exp(i7)(j(^)[jy], (25) 

and applying Eq. jSJ, it follows that 

2 

exp(i7) = ^4 — B ' (26) 

qABiN)q(^N)1iO)[N] 

and then a spin-weighted quantity V with spin-weight s transforms between the two patches as 

^fAT) = exp(is7)y(o). (27) 

There are many different cases for the transformation between the patches, and we give the transformation and 
Jacobian explicitly in only one case, when the inter-patch boundary is at 

TT TT 

P{0) ^, P{N) = Cr(o) = cr(Ar). (28) 

Then the coordinate transformation is 

f 1 \ / tan (7(0) \ 

Pm = arctan ( - ) , <J(n) = arctan ( j:^;^^ ) , (29) 

and the Jacobian evaluates to 



tanp(o)/ Vtanp(o) 



1^ = 1,^=0, (30) 



dcr^N) cosCTsinCT f^o'(7V) cos p sin p 



dp{o) cos^ (T sin^ p + cos^ p sin^ cr ' dcr^o) cos^ cr sin^ p + cos^ p sin^ cr 



(31) 



B. Computational implementation 

The existing stereographic code has been extended to the six patch coordinate system. Since the formulation of 
the equations in terms of spin- weighted quantities is independent of angular coordinate bases but dyad-dependent, 
it is necessary to re-implement only those objects that depend on the six patch dyad. We have therefore provided a 
numerical implementation of the new S-operators in the code. In addition, we have adapted the spin-transformation 
coefficients to the six patches. In order to test the six-patch code against the linearized solutions, the spin-weighted 
spherical harmonics sZ£m needed to be implemented for the six-patch coordinates and dyad. 

The S-operators are implemented via subroutines Di{s, e) and D2{s, ei, 62), which calculate 9 or B and the second 
derivative as the combinations 33, 33, 33 and 39. The parameter s specifies the spin s of the quantity to which the 
operator is applied, and ei, 62 G { — 1, 1} denote 9 and 9. 

The derivative operators d^, ffj. dndg, dp and da have been approximated by second- and fourth-order accurate 
centered finite difference stencils |2g. There is an input parameter that enables switching between second- and 
fourth-order accurate derivatives. 

In order to calculate the values of the 9-derivatives at the boundary of each patch, we need to access values from 
the neighboring patches. This is done by defining ghost zones, which contain the needed values of these patches. The 
multipatch infrastructure interpolates the spin-weighted (scalar) quantities between the patches using either cubic or 
quintic Lagrange polynomial interpolation, depending on whether angular derivatives are approximated by second- 
or fourth-order accurate finite differences. In both cases, the ghost-zone points lie on grid lines, so the interpolation 
was simple to implement because it was only one dimensional. After the interpolation, we apply the transformation 
law (|27() to transform quantities of spin-weight s 7^ to the current patch's coordinates and dyad. 

It turns out that we need a total number of 12 different spin-transformation coefhcients, since we have a P±- 
symmetry between the total number of 24 ghost zones across all Pj^i and P-i^ i = x,y,z patches. These coefficients 
are calculated and stored for repeated use in an initial routine. After each radial step during the integration of the 
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characteristic equations, the ghost zones are synchronized by the multipatch infrastructure, and afterwards, the code 
multiphes the appropriate spin-transformation coefficients with the synchronized values of the ghost zones. 

We have carried out checks of the angular grid including the Di- and _D2-operators and the spin-transformation 
coefficients, and we have found out that the code converges with second- and fourth-order accuracy, respectively. 

The spin-weighted spherical harmonics have been implemented by using the spherical harmonics in terms of the 
stereographic coordinate C = P + *9 ^i^d by applying the pseudo-numerical operators Di and D2 in order to obtain 
sZim for \s\ > 0. With pseudo-numerical, we mean that we apply the fourth-order Di and D2 operators with a very 
small delta-spacing such that we reach machine precision, since we know the Zfm everywhere and are not bound to 
the numerical grid. In order to use the stereographic routines for the Z^m, we have transformed the stereographic 
coordinate C to the six-patch coordinates and depending on the hemisphere, we use the stereographic routines for 
north or south patch, respectively. 

Furthermore, we have implemented an algorithm for calculating the linearized news ftmction fEg. ljlSI) ). in the form 

TV = lim ( fcll J _ 1 + d^p) . (32) 



IV. THE CONSTRAINT EQUATIONS 

If the boundary data satisfies the constraints i?oa = (here we restrict attention to the vacuum case), and provided 
the hypersurface and evolution equations are satisfied, the Bianchi identities guarantee that the constraints are satisfied 
throughout the domain Thus from an analytic viewpoint, evaluation of the constraints is redundant, but in a 

numerical simulation their evaluation may provide useful information concerning the reliability of the computation. 

We have written code that uses the Bondi-Sachs metric variables and derivatives to evaluate the following quantities 

i?oo, Rai and q'^RoA- (33) 

The expressions for the above quantities are very long and are not reproduced here. The Fortran code for these 
expressions was generated directly from the computer algebra (Maple) output. 

V. TESTING THE CODE 

In this section we first specify the linearized solutions against which the code will be tested, as well as the various 
parameters that describe a numerical solution and its output. Then we present, as Figures and Tables, the results of 
testing the comparative performance of the second order six-patch, fourth order six-patch and stereographic codes. 

We refer to the linearized solutions summarized in Sec. Ill CI In all cases we take 

v = 1 and m — 0. (34) 

We present results for the cases £ = 2 and £ — 3 with 

Ci = 3 • 10-^ C2 = 10-^ /3o = j ■ 10"^ {£ = 2) (35) 
Ci = 3 • 10-^ C2 = i • 10"^ /3o = i ■ 10"*' {£ = 3) (36) 

in Eq. 1)13(1 in the case £ — 2, and in Eq. p4|l in the case £ = 3. 

All the numerical simulations use a compactified radial coordinate x = r/{rwt +*") with r^t = 9. Data is prescribed 
at time it = as well as at the inner boundary r = 2 (which is equivalent to a; = 0.1888). The stereographic grids 
(with ghost zones excluded) are 

Coarse: = rig = rip = 41, Fine: rix = Uq = rip = 81; (37) 

and there is no overlap between the two patches, i.e. we set the code parameter ^gj^c ~ ^ which means that on the 
nominal grid, the holomorphic coordinate function = q + ip takes values in q,p € [—1, 1]. The six-patch grids are 
such that, over the whole sphere, the total number of angular cells is equivalent. We take 

Coarse: Ux = 41, n„ ~ Up ~ 24, Fine: Ux — 81, n„ — 7ip — 47. (38) 
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TABLE I: Comparative performance in the i = 2 case 



Stereograpliic six-patch, 2rid order six-patch, 4th order 

Averaged convergence rate of J 3.8456 3.8286 3.9112 

Averaged absolute error of J 3.3039x10"" 1.5491x10"'' 6.9157x10"" 

Averaged convergence rate of 3.3119 3.9642 3.5528 

Averaged absolute error of AT 2.2785x10"* 1.0913x10"* 8.4414x10"^° 

Averaged convergence rate of Tioo 1.2487 1.5000 2.0319 

Averaged absolute error of i?oo 3.1942x10"^^ 2.8779xl0~" 5.7668x10"^° 

Averaged convergence rate of i?oi 3.5560 3.5936 3.1296 

Averaged absolute error of i?oi 3.9214x10"" 1.6988x10"" 2.7331x10"^^ 

Averaged convergence rate of i?oA 3.4285 1.7558 2.0043 

Averaged absolute error of i?oA 5.2549x10"" 6.6397x10"" 2.1543x10"" 



TABLE II: Comparative performance in the 1 = ?, case 



Stereographic six-patch, 2nd order six-patch, 4th order 

Averaged convergence rate of J 3.9783 3.9106 4.0777 

Averaged absolute error of J 4.6461x10"" 3.2784x10"" 1.3677x10"^° 

Averaged convergence rate of 2.1201 3.9134 3.6262 

Averaged absolute error of A?" 4.9174x10"* 2.8182x10"* 1.7996x10"" 

Averaged convergence rate of iZoo 1.2743 1.7963 2.0330 

Averaged absolute error of i?oo 7.2594x10"" 4.5824x10"" 1.1744x10"" 

Averaged convergence rate of i?oi 3.5144 3.5383 3.3824 

Averaged absolute error of i?oi 1.3262x10"^" 7.5501x10"" 6.0924x10"^^ 

Averaged convergence rate of iioA 3.4326 1.9510 2.0156 

Averaged absolute error of i?oA 9.0299x10"" 1.0076x10"* 2.9654x10"" 



Six-patch results are reported for both second-order and fourth order differencing of the angular derivatives. In all 
cases, the fine grid has Au ~ 0.0125 and the coarse grid has Au = 0.025. Runs are performed for two complete 
periods, i.e. starting at u = and ending at w = 47r. 

Results are reported for the errors of the quantities shown using the L-i norm, evaluated at the time shown, averaged 
over all non-ghost grid-points over the whole sphere and between the inner boundary and future null infinity. The 
norm of the error in the news is averaged over the whole sphere at future null infinity. 

The error e for a quantity 4" is defined as 

*^(^) II ^n-fimeric ^ana/yiic ||- (^^) 

The convergence factor C is then defined as the ratio between the error e of low and high resolution 

m = (40) 



VI. CONCLUSION 



We have implemented a version of the characteristic numerical relativity code that coordinatizes the sphere by 
means of six angular patches. Further, the six-patch code has been implemented for both second-order and fourth- 
order accurate finite differencing of angular derivatives. 

We compared the errors in the second-order six-patch, fourth-order six-patch and stereographic versions of the 
code, using exact solutions of the linearized Einstein equations as a testbed. This was done for a variety of cases 
and using several different indicators to measure the error. The convergence rate of the metric (i.e. of J) was always 
approximately second-order, but in some cases we observed degradation of the order of convergence of other quantities. 
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FIG. 2: Convergence factor C and error norm e of J plotted against time in the case £ — 2 



all of which contain second derivatives of the metric. This has shown up in previous performed runs of the code and 
might be related to high-frequency error modes coming from angular patch interfaces or corners. On average, the 
error norm of second order six-patch was smaller than that of stereographic by a factor of order two (although there 
were cases in which the error was slightly larger). However, the fourth-order six-patch scheme exhibited a dramatic 
reduction in the error norm, by a factor of up to 47 compared to that of the stereographic case. 

Thus, we expect the six-patch characteristic code, in particular the version that uses fourth-order accurate angular 
finite differencing, to give significantly better performance than the stereographic version. 
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